#include <iostream>
#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp"

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 3>::tiny floatTinyArray_t;
typedef TunaArray<double, 3>::huge ScalarField3D;

int main () 
{

    double dt = 0.0001;
//
// Middle plane construction
//
    int Npoints_x = 5, Npoints_y = 5;
    floatTinyArray_t plane[Npoints_x][Npoints_y];
    double left_extreme = 0.1;
    double right_extreme = 0.9;

    double factor_x = (right_extreme - left_extreme) / (Npoints_x - 1);
    double factor_y = (right_extreme - left_extreme) / (Npoints_y - 1);

    cout << "\n ----- Plane construction ----- \n\n";
   
    for(int i = 0; i < Npoints_x; i++)
	for(int j = 0; j < Npoints_y; j++)
	{
	    plane[i][j](0) = left_extreme + i * factor_x;
	    plane[i][j](1) = left_extreme + j * factor_y;
	    plane[i][j](2) = 0.5;
	    cout << " plane[" << i << "," << j << "] = " 
		 << plane[i][j] << endl;
	}

    cout << "\n ----- Tracking: begins ... \n\n";

    StructuredMesh< Uniform<double, 3> > mesh(1, 47, 1, 47, 1, 47);
    floatTinyArray_t deltas = mesh.getDeltas();

// Velocity field at time t
    ScalarField3D u(mesh.getExtentNodes());       
    ScalarField3D v(mesh.getExtentNodes());
    ScalarField3D w(mesh.getExtentNodes());

// Velocity field at time t + dt
    ScalarField3D un(mesh.getExtentNodes());       
    ScalarField3D vn(mesh.getExtentNodes());
    ScalarField3D wn(mesh.getExtentNodes());

    ParticleTracking< RungeKutta4< Linear <double, 3> > > particle;

// Read the velocity at time t from a file
    InOut::readFromFile_2(u, v, w, 1584, "./Velocity/velc.");

// Read the velocity at time t + dt from a file
    InOut::readFromFile_2(un, vn, wn, 1585, "./Velocity/velc.");
//
// Now track every particle on the plane
    for(int i = 0; i < Npoints_x; i++)
	for(int j = 0; j < Npoints_y; j++)
	{
	    plane[i][j] = particle.track(plane[i][j], deltas, 
					 u, v, w, un, vn, wn, dt);
	    cout << " plane[" << i << "," << j << "] = " 
		 << plane[i][j] << endl;
	}

    cout << "\n ----- Tracking: ends ----- \n\n";

    return 0;
}



